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Abstract: We provide a bestiary of public codes and other algorithmic tools that can 
be used for analysing supersymmetric phenomenology. We also describe the organisation 
of the different tools and communication between them. Tools exist that calculate super- 
symmetric spectra and decay widths, simulate Monte Carlo events as well as those that 
make predictions of dark matter relic density or that predict precision electroweak or b- 
observables. Some global fitting tools for use in SUSY phenomenology are also presented. 
In each case, a description and a link to the relevant web-site is provided. It is hoped that 
this review could serve as an "entry-gate" and map for prospective users. 
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1. Introduction 

Analysis in high energy particle physics is becoming increasingly complex; the higher en- 
ergies and luminosities of current-day colliders lead to higher multiplicities in events. The 
current high-energy frontier is dominated by hadron-hadron colliders, at the Tevatron (pp 
at 2 TeV) or, in the near future, the Large Hadron Collider (pp at 10 or 14 TeV), lead- 
ing to additional complications in describing hadronic initial states and radiation. On the 
theoretical side, the currently most popular solution to the technical hierarchy problem is 
the Minimal Supersymmetric Standard Model (MSSM). A low energy parametrisation of 
the MSSM contains over 100 parameters. In fact, a truly supersymmetric version of the 
Standard Model contains one less parameter than the Standard Model, since the quartic 
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Figure 1: Schematic of the interaction between various programs that perform different SUSY 
phenomenology calculations. The need for information exchange is denoted by a line. 



Higgs coupling becomes a function of the electroweak gauge couplings in the supersym- 
metric version, instead of being a free parameter. However, in order for the MSSM to be 
phenomenologically viable, supersymmetry (SUSY) must be broken, and it is in the SUSY 
breaking sector of the model that the majority of the free parameters lie. The vast majority 
of this 100+ dimensional parameter space is ruled out by fairly tight constraints on flavour 
changing neutral currents. This is often taken to be evidence of some additional structure 
of the model in the flavour sector. High energy boundary conditions on the supersymmetry 
parameters that are flavour universal are popular, and may be motivated by various string 
(and/or grand- unified theory) models. 

A schematic of SUSY phenomenology calculations is shown in Fig. Typically, one 
may want to assume some high energy theoretical boundary condition upon the SUSY 
breaking sector. One wishes to calculate the MSSM spectrum and couplings consistent 
with this boundary condition and some input observables (Mz, mt . ■ .) with a spectrum 
calculator. The spectrum and couplings can then be passed to another program that 
calculates decays of the various sparticles. Once the masses and decays of the sparticles 
are calculated, this information may be passed to an event generator in order to randomly 
simulate several events in some high energy collision. This process is often split into two 
sub-steps: one performing the hard 2 — * N particle collision (matrix element generation), 
and one performing hadronic showering, initial state radiation and decays of the sparticles 
(event generation) . Experimental colleagues often then want to pass such simulated events 
through a detector simulation in order to see how the detector might smear the kinematics. 
Alternatively, the spectrum and coupling information could be passed to packages which 
calculate indirect observables, such as the dark matter relic density left in the universe, 
dark matter direct detection cross-sections, electroweak or 6-observables. These data may 
then be used in global fits to the particular SUSY breaking scenario assumed. Some of 
the programs available perform several of these tasks, but there is currently no single 
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program that performs all of the tasks. Previously, information was passed around on an 
ad hoc basis: each spectrum generator had to be interfaced separately with each program 
that used its output. With N independent codes, the required number of interfaces such 
that they could each exchange information to the others was ~ Nl. For this reason, 
several accords have been written and agreed upon in order to cut down on the total 
number of required interfaces, with an associated reduction in the number of mistakes in 
the interfacing procedure. 

The SUSY Les Houches Accord (SLHA) [|IJ allows information on the masses and de- 
cays of SUSY (and some relevant Standard Model) particles to be passed in between codes. 
The accord is based on ASCII text, in order to allow easy cross-language communication 
without introducing platform dependence. The parsing of (files or memory variables con- 
taining) such ASCII text is an easy task for many human beings, but the disadvantage 
of an ASCII format is that developers of tools must write parsing code. Luckily, even 
this task has been performed, with a SLHA-file parser available Pj. The original SLHA 
dealt purely with the "vanilla-MSSM" : inter-generational sparticle mixing is not taken into 
account, R-parity and CP are conserved. The second SLHA ||] generalises the possible 
MSSM models: R-parity violating, CP and flavour-violating versions of the MSSM are 
all specified. In addition, the most popular MSSM extension where a Standard Model 
singlet chiral superfield is added, the Next-to-Minimal Supersymmetric Standard Model 
(NMSSM), is covered. 

The original Les Houches Accord (LHA) Q, allows hard parton- level events to be 
passed from matrix element calculators onward down the chain to the event generators. 
It does this by means of a fortran77 common block, which specifies properties of the 
particular process being simulated such as the types of particles involved and their mo- 
menta. Colour flow in the diagrams requires particular attention and is encoded in the 
LHA. However, all of the Les Houches Accords attempt to hide such details and require- 
ments from the user. Only tool developers have to concern themselves with them. More 
recently, the Les Houches Accord event record has been changed to a minimal XML-style 
structure, for clarity, simpler parsing and to side-step cross-language difficulties || and 
several parsers (in different languages) have been developed and are available. The accord 
has also been re-written to take into account potential new beyond the Standard Model 
physics models ||. 

In this review, I shall briefly describe the publicly available, supported, documented 
codes which allow supersymmetric phenomenological calculations 1 . In each case, a link to 
a current web-site and a reference to the relevant manual is given. The default language 
of each program is f ortran77, but if the code is written in a different language, it shall be 
detailed in this review at the point when the main functionality of the code is discussed. 
As time passes, it is foreseen that some of the links listed here will become out of date. The 
reader is advised to read the manual of any code they wish to use from the electronic arXiv 
web-site in order to find updated links to downloads etc. In addition, more accuracy and 
extended functionality will no doubt be added to the various programs as time passes. This 

1 In fact, a "quick guide" of SUSY tools was written over two years ago The present review contains 
an updated and much more extensive overview of the field. 
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guide is intended to serve as a snap-shot of documented, supported, publicly available SUSY 
phenomenological tools at the time of writing. It is not practical to continually up-date it 
as the state-of-the-art evolves. However, it should be mostly accurate for a good few years 
and there are plans to extend a Beyond the Standard Model tools repository ||] to a Wiki 
site, so that the authors of the codes may include up-dates to the accuracy/functionality 
as they occur. We shall not describe here any of the detector simulations. First, in section [2], 
we shall describe codes that calculate MSSM SUSY spectra and decays. Then, in section |3|, 
we list matrix element generators, followed by event generators in section ||. We then turn 
to constraints: in section [| we review public SUSY dark matter codes, followed by other 
indirect constraint calculators in section ||. We review some of the algorithms required 
to perform global fits to SUSY models using indirect observables in section 0. Finally, in 
section |8[ we conclude and present a table summarising the functionality of the SUSY tools 
mentioned in this review. 

2. Spectrum and Decays 

There are four publicly available dedicated MSSM spectrum generating codes, displayed 
in Table |l|. They all solve the MSSM renormalisation group equations (RGEs) to two-loop 
order, subject to two sets of boundary constraints. One set of boundary constraint is at the 
weak scale, and matches the MSSM parameters to current data on Standard Model particle 
masses and couplings. It also ensures that electroweak symmetry is broken successfully 
by adjusting the MSSM fi parameter. The other boundary condition is typically at a 
high energy scale, and involves setting the SUSY breaking parameters according to some 
theoretical model of SUSY breaking mediation. Universal mSUGRA, minimal anomaly 
mediation and minimal gauge mediation are supported by all of the codes. In addition, 
non-universal models such as those that can be invoked by the SLHA are supported. Each 
of the codes supports different additional possible SUSY breaking models. They each also 
support the SLHA aside from ISA JET ||. An unofficial version of ISAJET which outputs 
in SLHA format does exist, however |fic| l. 



Name 


Language 


RGEs 


comment 


manual 


ISAJET 
SOFTSUSY 

SPheno 
SUSPECT 


C++ 
f ortran90 


dominant 3rd 
3-family mixing 
3-family no-mixing 
dominant 3rd 


VR 
VR 


I 
I 


B 

ii 

12 
13 


] 
] 
1 



Table 1: SUSY Spectrum generators, vr indicates that the program includes an option for in- 
cluding right-handed neutrinos in the spectrum in order to obtain neutrino masses. Dominant 3rd 
RGEs mean that all Yukawa couplings other than the third family's are neglected in the RGEs 
whereas 3-family no-mixing means that all diagonal Yukawa couplings are included. 

The details of approximations within the codes can be found within the manuals, and 
although similar, do tend to vary somewhat. They may differ by higher-order correc- 
tions, for example. The matching conditions to current data at the weak scale is mostly 
in the one-loop approximation. But when one is correcting a QCD cross section with 
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Comparison of SUSY spectrum generators: 
mass spectra, relic densities, etc 

Updated on 27 Feb 06 to Isajet7.74 and SoftsusyZ.O.S 

On this website you can compare the mass spectra of four public 
SUSY codes, sajet7.74, Softsusy2.0.5 J Spheno2.2.3 and Suspect2.3.4, 
together with the resulting relic densities as computed by micrOMEGAsl.3 

Please note: 

• Isajet has mb(mb)=4.214 alpha s(MZ)=0.1172 hard-coded in the program 

so any change of mb(mb) and aTpha s below will MOT apply to Isajet 
» The improvements in Isajet7.74 are based on hep-ph/0511123. 

NEW: you can now also download the corresponding SLHA spectrum files! 

mSUGRA input: m_0 [Too GeV 

m_l/2 |2M GeV 

A_0 |Too GeV 

tan beta [To (ca 1.6 -50) 

sign(mu) |+1 

SM input: m_t [L75 GeV (anshell value) 

m.b(mb) [OTS GeV (MSbar) 

alpha_s( MZ) |0.1172 ( MSbar) 

submit | reset | 

Figure 2: Comparison web-page mSUGRA form p3 



stop loops, for example, in order to extract the MSSM value of a s (Mz) that should be 
used, the question arises which stop mass should be used? The codes either use pole 
masses or running masses evaluated at different scales for the stops running in the loops. 
The difference between these choices is actually a higher order effect, and could only be 
fixed were the two-loop matching conditions known. Broadly speaking, RGE evolution is 
two- loop, in different approximations, as displayed in Table [l]. ISAJET differs from the 
other codes in that it decouples sparticles at their mass scales, thus re-summing terms 
~ C(l/(16vT 2 ) log[AM/M]), where AM is the splitting between two sparticles and M is 
their average mass. On the other hand, the other three codes all evolve using MSSM RGEs 
above Mz without decoupling sparticles, but then one-loop decoupling effects are added 
to the weak-scale boundary condition to leading logarithmic order. This latter approach 
allows the easy addition of some one-loop finite pieces, some of which are missed by the 
mass-scale decoupling approach taken by ISAJET. In summary, one may expect the mass- 
scale decoupling approach to provide a more accurate answer when sparticle splittings are 
very large, and Mz decoupling including all finite terms to be more accurate for a more 
typical sparticle splitting. The codes all agree to the percent level, except in particularly 
difficult parts of parameter space such as the focus point or very large tan/3 [fig ], where the 
differences can be much larger due to inherent numerical instabilities and the size of higher 
order corrections in those regions. Fortunately, a web-site exists |]T3 where one can input 
a SUSY breaking point on a web form as exemplified in Fig. [2], and quickly compare the 
output from the different codes. If one is doing a study on a particular point, for example, 
this provides a quick practical way of finding out if the point comes with particularly large 
theoretical uncertainties or not. 
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2.1 MSSM Sparticle Decays 



Currently, the programs that calculate sparticle decay branching ratios are Herwig++ [16] 



(C++), ISA JET @, MadGraph 0, PYTHIA |]|, SDECAY p| and SPheno [|T|]. SDECAY, 
PYTHIA and Herwig++ take the SLHA stream from any of the other codes in order to pro- 
duce SLHA-compliant output including decay information, whereas the other two codes are 
linked to their spectrum generators. The decay packages implement tree-level two-body de- 
cays of fermions and gauginos and three-body decays of charginos, neutralinos and gluinos. 
SPheno includes gluonic QCD corrections into decays by quarks. SDECAY implements some 
three and four-body decays of top squarks and one-loop corrections to the two-body decays. 
PYTHIA and Herwig++ contain internal routines for calculating sparticle decays, including 
tri-linear R— parity breaking effects. Herwig++ and MadGraph include angular correlations 
between subsequent decays in a sparticle cascade decay using the pioneering techniques 
of ref. |20|| , whereas all the other codes make a phase-space approximation. The program 
BRIDGE [^] was written in order to decay particles passed to it by matrix element gener- 
ators in general models defined in the MadGraph format, then pass them on to showering 
and hadronisation programs. It calculates two and three-body tree- level decays itself, while 
keeping track of initial vertex spin structures via HELAS. Typically phase-space is a reason- 
able approximation in hadronic collisions unless one is trying to fit the spin of sparticles. 

2.2 Higgs Masses and Decays 

There are some packages specialising in SUSY Higgs calculations: FEYNHIGGS |22j calculates 
the Higgs masses in a Feynman diagrammatic approach. In calculating Higgs masses, 
important two-loop effects are included for the MSSM with or without complex parameters, 
with a re-summation of the leading (s)bottom corrections. One-loop non-minimal flavour 
violating corrections to Higgs masses/mixings are included at the one- loop level. The 
program calculates the Higgs spectrum and decays and provides an estimate of theoretical 
uncertainties in the prediction. The two-body tree level decays include dominant one-loop 
corrections and the Higgs decays to gg and 77 include all of the MSSM particles in the loop. 
A FEYNHIGGS web-form interface exists for checking single points in parameter space p2|. 



The program CPsuperH [ 23 1 also performs MSSM Higgs calculations when CP violating 
phases are present, including some effects up to two-loop order. The program is based 
on renormalisation-group improved diagrammatic calculations that include logarithmic as 
well as threshold corrections and 6-quark Yukawa coupling re-summation. Some dominant 
one-loop pieces are included in the Higgs decays, which can be into SUSY or SM particles 
(including some important three body decays). The Higgs couplings and neutral Higgs 
mixings are also provided by the program. HDECAY 1 24] calculates up to three-loop QCD 



corrected decays of Higgs bosons in the CP-conserving MSSM where expressions exist in 
the literature, including some loop-induced decays, decays into two massive gauge bosons, 
three-body decays and decays into SUSY particles. Leading electroweak corrections are 
included (they can become important in the large Higgs mass regime due to enhanced Higgs 
self- interactions). All MSSM particles are included in the loop for the calculation of 77 and 
gg Higgs decay modes. The leading QCD corrections are included for the gluonic mode. 
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FCHDECAY [25j computes the flavour changing neutral current decays BR(H° — * tc,tc) 
and BR(H° -> 6s, 6a) in the flavour violating MSSM, using SLHA2 for input/output. It 
includes full one-loop SUSY QCD contributions. 

2.3 NMSSM 

The addition of a Standard Model singlet superfield to the MSSM constitutes a potential 
solution to the \i problem (why fi is of order the electroweak scale rather than some much 
heavier scale) and is called the next-to-MSSM (NMSSM). In the package NMSSMtools 
sparticle masses are calculated using two-loop NMSSM RGEs in the dominant third family 
approximation. Tree level sparticle decay widths and branching ratios are also calculated. 
The Higgs masses, couplings and widths (for two-body modes) are calculated within the 
NMSSM using approximations to the one and two-loop dominant corrections. For decays 
into the SM particles, the widths are calculated including one-loop SM QCD corrections. 

3. Matrix Element Generators and Cross Section Calculators 

In the high-energy LHC regime, often we wish to calculate the production of more than two 
hard particles. This is the job of matrix element generators. Matrix element generators can 
usually calculate total or differential cross-sections and/or produce independently sampled 
events. Simulating (for example) the production of two squarks plus some additional hard 
QCD radiated jets requires us to deal with complicated Feynman diagrams involving many 
particles in the final state. For this job, one uses a matrix element generator, which simu- 
lates or calculates the hard process (e.g. qq — > t\t\ + jet). The matrix element generators 
are currently mostly at tree-level, particularly as regards SUSY physics. In practice, 2^6 
to 2 — > 8 processes may be feasible depending upon the number of Feynman diagrams, 
although a vast amount of CPU time may be needed to compute them (using, for example, 
the grid). The number of Feynman diagrams tends to grow to be too large with increasing 
numbers of final-state particles. 

FeynArts/FormCalc are Mathematica packages for the generation and calculation 
of Feynman diagrams up to one-loop order. They can thus be used to calculate matrix 
elements for scattering processes. Up to 2 — > 3 processes can be calculated at the one- 
loop level with integration optimisation, although FormCalc has been successfully used to 
compute 2 — * 4 processes at tree-level. Vanilla, CP violating and non-minimal flavour 
violating versions of the MSSM are available. There is also a way of encoding some new 
physics Lagrangian model for extensions. FormCalc simplifies the amplitudes generated 
by FeynArts analytically and generates fortran77 code for the numerical evaluation of 
the squared matrix element. Automatic generation and pictorial representation of Feyn- 
man diagrams is also supported, as is convolution with parton density functions (PDFs). 
Recently, a program HadCalc |28| has been developed based on FeynArts and FormCalc. 
It takes the output from those codes in terms of partonic cross sections and convolutes 
them with PDFs. There are convenient ways to place cuts and an interactive menu-driven 
front-end that can be used to dial in SUSY parameters. 
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Currently, CalcHEP [29| and CompHEP [30] (C) 2 can cope with up to 6 external legs in 
a Feynman diagram, for example (1 — > 5 or 2 — > 4). The two programs can produce C 
output of analytical expressions for subsequent compilation and use. They have graphi- 
cal interfaces which can display modulus squared Feynman diagrams. Models are already 
defined for the MSSM, NMSSM and the CP-violating MSSM. For the encoding of model 
Lagrangians and parameters, LanHEP [31] (C) is used. If the user wishes to extend some 
SUSY model outside of the ones already defined, LanHEP provides the means. MadGraph 
performs vanilla MSSM matrix element calculation with SLHA input. Helicity amplitudes 
are constructed based on the HELAS |3^] library in order to encode spin information of the 
produced particles, which can be used in their decay. Feynman diagrams are drawn, and 
fortran77 output is produced for the matrix element. A Monte-Carlo integrator pack- 
age has been included in MadGraph, and SUSY differential or total cross sections can be 
calculated using it. Alternatively, the final result can be Les Houches Accord formatted 
parton-level events that can be passed into an event generator for subsequent parton show- 
ering and hadronisation. The MadGraph web-site has a form that can be filled in to get 
events returned automatically. The event generator SHERPA |33| utilises an event generator 
Amegic++ |34| (C++) that also uses the helicity amplitude technique and calculates at the 
tree-level, with the possibility of up to six particles in the final state of the hard scattering. 
Whizard [35 (f ortran95) includes initial and final state polarisations and can calculate in 
the vanilla MSSM as well as the CP-violating case. It uses O'Mega [36| (O'Caml) to trans- 
late a helicity amplitude into computer code as needed. O'Mega is designed with special 
tricks to avoid the factorial increase in CPU time with the number of external particles. It 
has been demonstrated to work for some processes with eight particles in the final state. 
SUSYGEN [[}?]] is restricted to 2 — > 2 SUSY production processes. It can include polarisation 
in e + e~ collisions and covers vanilla as well as R-parity or CP-violating MSSM models. 
GRACE |3£| performs computations of e + e _ — > up to four bodies in the MSSM at tree-level. 
GRACE draws the relevant Feynman diagrams for the user. 

Numerical results of several hundred SUSY production cross-sections were compared 
between MadGraph 0, SHERPA Q and Whizard Q and they were all found (eventually) 
to agree @. PR0SPIN0 @ (fortran90) computes MSSM next-to-leading order cross 
sections for the production of two sparticles at hadron colliders. It can also cope with the 
production of weak gauginos in the split SUSY framework. Detailed calculations of cross- 



sections of e + e — > sleptons at the one-loop level are also available from ILCslepton [41 



FeynHiggs [22] calculates Higgs production cross-sections for the Tevatron and the LHC 
including SUSY corrections at the production vertex. 



4. Event Generators 

The most well-known general purpose SUSY event generators (PYTHIA jl^], Herwig++ fill ) 
usually implement a hard-sub process in terms of two particles scattering on two particles, 
represented by the central vertex in Fig. ||. The initial particles in this hard-sub process 

2 These two programs have the same origin, but at some stage the development of them branched. 
Because of this, although the programs are now different, many features of them are similar. 
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Figure 3: Schematic of a hadron collision simulation in an event generator. 



may be leptons, or point-like constituents of hadrons. In the latter case, the quarks and 
gluons are extracted from a hadron by means of the parton density functions. The hard- 
sub process is usually calculated at leading order in perturbation theory within event 
generators, especially for exotic signals such as SUSY. If SUSY particles are produced in 
the simulated event, they are then decayed randomly, according to the branching ratios 



calculated by a program in section 2A. The resulting cascade decay spits out SM particles, 
some of which may be quarks or gluons. These quarks and gluons then emit soft QCD 
radiation, which is modelled by the parton shower. Parton showering encodes the fact that 
the matrix elements of massless coloured particles emitting a gluon have a singularity in 
the infra-red or collinear limit. The initial state may also shower, emitting QCD radiation. 
It can be important to include effects in the shower coming from colour coherence in order 
to describe the resulting jets adequately. Various properties such as angular ordering of 
the shower (with preceding emissions being at smaller angles) are evident in the resulting 
event. Once the partons are showered down to some energy scale to be decided by the 
event generator, some non-perturbative modelling (and a tune to data) collects the partons 
together into hadrons which, after their decays have been simulated, may be observed in 
the detector. Finally, the simulated events are often represented by a series of lines in a 
text file (representing variables held in memory), each line describing the kinematics and 
state of a particle involved with the event. Information on which particles decayed into 
which other particles is also indicated in this event record. We briefly mention some of the 
Standard Model properties of event generators, since elements of them are also relevant to 



SUSY events, but the interested reader should see ref. [42] for a more complete guide to 
Standard Model event generators in hadron collisions. 

A general framework for encoding new physics models is included in Herwig++ jjlfj] so 
that users may define the relevant particles and Feynman rules for the hard sub-process. 
The MSSM has already been defined within the framework, but extended SUSY models must 
be input by using it. It can read SLHA files for input SUSY information. The Herwig++ 
shower algorithm treats QCD radiation from coloured heavy objects (for example tops). 
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It can evolve to zero transverse momentum of emissions, giving an improved simulation of 
the dead-cone effect for radiation from massive particles. An eikonal multiple scattering 
model is used for simulating additional partonic collisions in the same hadronic collision. 
Such processes form part of the underlying event. Herwig++ uses a cluster model for the 
hadronisation step, clustering quarks and gluons that have similar kinematics into colour 
singlet states, which decay to hadrons and hadron resonances. Herwig++ treats its decays 
including spin correlations all of the way down the various decay chains. Herwig++ relies 
on an underlying C++ structure developed for high energy collisions called ThePEG [43]. 

PYTHIA [18], on the other hand, does not depend upon ThePEG, but is stand-alone. 
As well as hadron-hadron collisions, it can deal with e + e~ beams. Initial and final-state 
parton showers are based on pr-ordered evolution, terminating at 1 GeV. Although there 
is currently a C++ version in development, it does not contain any SUSY physics and so we 
concentrate on the older f ortran77 version in this review. Many different options for the 
changing the models of parts of the PYTHIA simulation are possible, but here we describe 
the default models. Hadronisation and hadron fragmentation (decay) are modelled by the 
Lund string model, where hadrons are modelled to be a colour flux tube, ended where 
the (di-)quarks are located. The MSSM, the NMSSM, tri-linear R— parity violation, as 
well as long-lived coloured sparticles such as those that exist in models of split SUSY, are 
included in the PYTHIA distribution. Polarisation is included for e + e~ incoming beams. 
All decays of sparticles are using the phase-space approximation, and so sparticle spin is 
not simulated. 

A new event generator has recently been developed called SHERPA [|33| (C++). The 
main design feature of SHERPA is that it combines parton shower evolution and matrix 
element generation. The MSSM with or without CP violation and full inter-generational 
mixing is included, as well as a general formalism (compatible with the matrix element 
generator Amegic++ [I34], which ships with the SHERPA distribution) provided for adding 
new particles and interactions. Amegic++ can practically handle up to six jets in the 
final state for e + e~ collisions, and up to three jets for hadron collisions. One of the 
difficulties of combining parton showers and matrix element generation for hard jets is 
the problem of double-counting. If one simply adds the matrix element generation for 
3 jets to the 2-jet plus parton shower sample, one could easily double count the region 
where one of the jets is soft (and therefore already included in the parton shower). In 
SHERPA, the parton shower evolution and matrix element generation are matched via the 
CKKW formalism [44], where the matrix element configurations are re- weighted according 
to a pseudo shower history and shower emissions that overlap with higher order matrix 
elements are rejected. SHERPA also performs the hadronisation/fragmentation step using 
its own cluster model [45], which includes di-quark spin effects and a dynamic separation 
of the regimes of clusters and hadrons according to their masses and flavours. When 
calculating SUSY decays in SHERPA, there is currently no facility for picking the decay 
products automatically, the user must supply which decay chain is required. After this has 
been done though, the spin information and off-shell effects are included in each sparticles' 
decay into the next sparticle and Standard Model particle. 

The ISAJET B event generator simulates pp, pp and e + e~ collisions at high energy, 
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based on perturbative QCD and phenomenological models for parton and beam jet frag- 
mentation, not including colour coherence effects: the probability of emitting a soft gluon 
is multiplied by a factor given by the Alterelli-Parisi function. ISAJET keeps only the parts 
that are in the exact collinear limit, but uses non-collinear kinematics. QCD radiation 
from initial and final states is simulated. Sparticle pair production at tree-level is sup- 
ported, along with subsequent decay. The ISAJET hadronisation model is the independent 
fragmentation ansatz of Field and Feynman, which forms new (di-) quark- anti-quark pairs 
out of partons, and groups them together into mesons and baryons with some fraction of 
their summed momenta. 

5. Predictions of SUSY Dark Matter 

The recent WMAP5 cosmological fits to the cosmic microwave background and other data 
provide us with an accurate observation of the density of dark matter in the universe as a 
fraction of the relic density: Qh 2 = 0.1143 ± 0034 @. The MSSM offers several possible 
candidates for weakly interacting massive particles that could play the role of cold dark 
matter, since the lightest supersymmetric particle (LSP) is stable from the assumption 
of R— parity. The dark matter candidate obviously must be without electric charge, so 
that it does not interact with light, and also should be colourless, otherwise it would have 
fused with nuclei during nucleosynthesis and been discovered in anomalously heavy isotope 
searches. Gravitinos and lightest neutralinos are possible candidates within the MSSM, 
although in extended models, other SUSY particles are possible dark matter candidates. 
In principle, SUSY dark matter may be discovered in direct detection experiments, where 
nuclear recoils from collisions with SUSY dark matter are possible. If the dark matter can- 
didate interacts too weakly (for example in the case of the gravitino) , the direct detection 
cross-sections are far too small to small to be seen in the foreseeable future. However, in 
the case of the neutralino, there is a chance for direct dark matter observation. In order to 
detect dark matter on earth, it is of course a necessary condition that there is some dark 
matter going through the detector. While only a small amount is empirically known about 
the small-scale structure of dark matter halos, numerical iV-body simulations indicate that 
even if one starts with strict filaments or cusps of dark matter, subsequent Newtonian 
evolution will tend to smear it out. Thus, the prospect of having our galaxy's dark matter 
localised completely elsewhere in the galaxy seems unlikely, but it should be borne in mind 
that any calculation of the local dark matter flux on earth is subject to large astrophysical 
modelling uncertainties. Aside from the direct detection of particulate dark matter, there 
are prospects for indirect detection, where for example, dark matter annihilation in the 
sun, in the earth or in the centre of the galaxy produces high energy particles that can be 
detected on earth or on satellites. 

Once a SUSY model's parameters has been fixed and a cosmological model is assumed, 
it is possible to estimate the amount of current dark matter relic density in the universe 
is predicted. The codes tend to assume the A CDM model of a cosmological constant plus 
cold dark matter component, since this is quite a simple and good fit to the WMAP and 
large scale structure data. One has to track the abundances of the different species of 
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sparticle through the early evolution of the universe. They can annihilate with each other 
into Standard Model particles, but eventually, the expansion of the universe makes the 
sparticles too far apart to interact. Aside from losses due to annihilation processes, each 
sparticle will end up as a LSP through its decays. The tracking of the abundances of the 
various sparticle species and involves the evolution of coupled Boltzmann equations. There 
are many different annihilation cross-section processes to consider, and the relevant public 
codes currently calculate at tree-level. The velocity distributions of the SUSY particles are 
derived from Maxwell-Boltzmann approximations. This can still involve the calculation of 
thousands of Feynman diagrams, however. As a consequence, the current tools calculate 
mainly at tree-level. However, loop corrections can give large ~ O(10)% effects in some 
cases [47|. 



DarkSUSY [48] contains hard-coded matrix elements for the many different annihilation 
processes of the vanilla MSSM. It can calculate the relic density as well as direct and 
indirect detection rates, with a choice of different nuclear form factors for the direct rates. 
Solar system WIMP velocity distributions can be used to calculate the capture in the 
Earth of dark matter particles. An exotic component of positron, anti-proton and anti- 
deuteron in cosmic rays originating from neutralino pair annihilation in the galactic halo 
can be calculated. Hadronisation and fragmentation was calculated with PYTHIA and the 
results tabulated from various neutralino masses, which DarkSUSY interpolates in order to 
provide an estimate of the particle yield. For particle yields coming from annihilation in 
the earth and the sun, 6 fundamental channels are included: cc,tt,bb,T + r~ ,W + W~ and 
Z°Z°. Recent solar and terrestrial density models are included as a necessary ingredient in 
the calculation. For galactic halo annihilations, W + W~, Z°Z°, W + H~, Z°h°, Z°H°, h°A° 
and H°A° channels are included, with subsequent decay of these particles, including the 
heavy quarks c, b and t. The gg, 77 and Z7 channels occurring at the one loop level are 
also included. Anti-matter production yields from dark matter annihilation in the Galactic 
halo are determined by DarkSUSY. Modelling the propagation of anti-matter is non-trivial, 
but DarkSUSY attempts this through various approximations which can be found in the 
manual. High energy neutrinos and neutrino-induced muons can be detected by neutrino 
telescopes and their yields are calculated. The SLHA is currently not supported, but 
instead dedicated interfaces to ISAJET and SUSPECT are included for spectrum generation. 
There is a web interface linked from the DarkSUSY homepage for inputting a MSSM model 
and calculating the relic density and detection cross-sections. Thus, if one is doing an 
analysis on one point in parameter space, one can check its dark matter properties easily 
on-line. 

IsaRED is part of the ISAJET || package, and calculates the relic density of neutralino 
dark matter in the MSSM. Annihilations between x?> X2> Xi > ^1, fix, f\, D e , D^, u T , u\, c%, 
t\, d%, s%, bi and gluinos are taken into account in the calculation. In the same package, 
IsaRES evaluates spin- independent and spin-dependent direct detection rates. Squark, Z° 
and Higgs exchanges are included at tree-level and neutralino-gluon interactions involving 
quarks, squarks and Higgs bosons are included at the one-loop level. 

micrOMEGAs (C) |^9| calculates the relic density of the LSP at tree- level and direct /indirect 
detection rates in the vanilla MSSM, the MSSM with complex phases and the NMSSM. 
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Important higher-order QCD and SUSY QCD corrections to Higgs quark vertices are in- 
cluded. The program can be used to calculate the relic density of a charged and/or coloured 
next-to-lightest supersymmetric particle. This can be useful in the case of a gravitino LSP. 
Gravitinos are not simulated by micrOMEGAs, but a simple formula can be used to extract 
their relic density from the relic density of the next-to-lightest supersymmetric particle. The 
most important annihilation channels for any given model point can be output. micrOMEGAs 
uses CalcHEP [29] in order to calculate any necessary Feynman diagrams, and so extensions 
can be encoded using the LanHEP [[H]] Lagrangian formulation, for models where there is 
only one stable particle. Only diagrams that may contribute up to some specified fraction 
(by default, 10 -5 ) of the thermally averaged total annihilation cross-section are included, 
which makes for faster computation. LSP scattering rates on nucleons and nuclei in the 
spin-independent and spin-dependent interaction cases are also presented. 7, e + , p and v 
yields for indirect direction purposes (at v — > and/or in the continuum) are calculated. 
Like DarkSUSY, micrOMEGAs uses the basic channels qq, t + t~, W + W~ and Z°Z° 

and interpolates tables of 7, e + ,p and v production as obtained by PYTHIA. For channels 
that contain two different particle species AB, the final spectrum is obtained by taking the 
average of di-A production and di-B production as a rough approximation. The galactic 
gamma-ray flux is calculated with a modified isothermal distribution of dark matter in 
the galaxy. For direct detection rates, higher order corrections to Higgs-quark vertices and 
one-loop neutralino-gluon interactions are included for the vanilla MSSM, the CP-violating 
MSSM and the NMSSM. 



6. Predictions for constraints 



There are many constraints upon supersymmetric models: direct constraints tend to be the 
easiest to implement, being (usually) phrased as lower bounds on sparticle masses. Rele- 
vant indirect constraints are upon branching ratios for rare decays, precision electroweak 
observables or electric dipole moments for example, and often occur at the loop level. Be- 
ing of general utility, FormCalc [^7]] can be used to calculate the relevant SUSY matrix 
elements. 



6.1 b observables 

The branching ratio of b — ► 57 has long been used to constrain supersymmetric models, and 
is calculated by several codes. Many of the codes calculate it in the vanilla MSSM without 



SUSY flavour mixing. For minimal flavour violating MSSM computations, SusyBSG [50] 
calculates the branching ratio for the decay b —* sj taking into account all of the avail- 
able next-to-leading order (NLO) contributions, including the complete supersymmetric 
two-loop QCD corrections to the Wilson coefficients of the magnetic and chromo-magnetic 
operators, as well as an improved NLO determination of the relation between the Wil- 
son coefficients and the branching ratio. micrOMEGAs f49[| , predicts the branching ratio 
including next-to-leading order contributions for the Standard Model. The charged Higgs 



and supersymmetric large tan/3 effects beyond leading-order are included. DarkSUSY [48] 



performs a NLO calculation which is complete for the Standard Model prediction and 
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adds some dominant NLO MSSM corrections. SPheno |Tj and SUSPECT |lj include one- 
loop MSSM corrections and some NLO corrections to the branching ratio. Superlso [51] 
calculates the b — > S7 branching ratio in the vanilla MSSM with flavor violation, NLO 
supersymmetric contributions and next-to-next-to-leading order (NNLO) Standard Model 
contributions. Flavour violation is supported through the SLHA2 interface. Superlso is 
currently the only code to predict the isospin symmetry breaking Ao_ of the B — ► K*"f 
decay including the NLO SUSY contributions. CPsuperH [23] can provide a prediction for 
the branching ratio as well as as its CP-asymmetry and SUSY contributions to B® d — B® d 
mass differences (AM S ^). FeynHiggs [^] provides a prediction for the b — > S7 branch- 
ing ratio including non-minimal flavour violating effects. Another 6— physics observable 
that can constrain SUSY is the to-date unobserved rare decay mode B s — ► /i + pT . The 



SUSY calculation in mi cr OMEGAS [49] includes the one-loop contributions due to chargino, 
sneutrino, stop and Higgs exchange, re-summation effects at high tan (3 are taken into 
account. CPsuperH also performs the calculation of BR(B S — > p + p~) in the CP vio- 
lating MSSM, as well as B& — > t + t~ , B u — ► t + v t . Each branching ratios is calculated in 
the single-Higgs insertion approximation. NMSSMtools |2g] calculates b — > S7, B s — > 
and B + — > t + u t branching ratios as well as AM S ^ in the NMSSM at one-loop order. 
ISATOOLS H includes NLO contributions to some of the Standard Model Wilson coeffi- 
cients for BR{b — > S7) and one-loop MSSM corrections. Branching ratios for B s — > 
and Bd — > t + t~ are calculated to one-loop, using approximations for the chargino masses 
(neglecting their mixing). The fitting program SuperBayes [^] uses the micrOMEGAs pre- 
diction of BR(b 57) at NLO and then augments it by NNLO Standard Model QCD 
contributions. 

6.2 Anomalous magnetic moment of the muon 

The anomalous magnetic moment of the muon is currently around 3<7 higher than the Stan- 
dard Model prediction. There is thus room for a non-zero SUSY contribution. ISATOOLS H] , 



SPheno |0, Superlso pi, micrOMEGAS pi and DarkSUSY [p| calculate the predicted 



SUSY contribution to one-loop order, whereas FEYNHIGGS |22] and SUSPECT [T^] also in- 
clude some two-loop corrections. 

6.3 Electric dipole moments 



For calculations of electric dipole moments in the CP-violating MSSM, micrOMEGAS [49| 
can provide estimates for the electron and Thalluim. One-loop neutralino/chargino con- 
tributions and two-loop squark, quark and chargino contributions are included as well as 
four-fermion operators for Thallium. Two-loop Higgs-mediated contributions to electron, 
muon and Thallium electric dipole moments are calculated in CPsuperH [23]. However, cur- 



rently some well-known one-loop contributions have yet to be implemented. The Thallium, 



neutron and mercury electric dipole moments are calculated in FEYNHIGGS [22] 
6.4 Electroweak observables 



micrOMEGAs |4!| and SUSPECT [ 13 1 can output the Ap parameter, which describes some 



loop corrections to electroweak observables. They both contain one-loop stop/sbottom 
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contributions, as well as two-loop QCD corrections due to gluon exchange and the heavy- 
gluino limit of gluino exchange. FeynHiggs [22 also contains a calculation of Ap, with 
corrections up to two-loops. SPheno [12] outputs the one-loop sfermion contributions to 
Ap. In terms of the electroweak observables themselves, FeynHiggs also computes M\y 
and sin 2 9w ^ including some two- loop SUSY contributions, non-minimal flavour violating 
effects and the effect of complex phases in the stop/sbottom sector at one- loop. 



7. Fitting tools 

We first introduce some necessary statistical terms, then go on to discuss their use in the 
context of SUSY fits. Typically, global fits of models to data utilise a statistical "figure 
of merit" for each point in parameter space to characterise how well it fits data. The 
most familiar one for particle physicists is probably x 2 > but sometimes likelihood is used 
instead. Likelihood L can be simply related to the x 2 parameter, L oc e~* 2 / 2 . L or x 2 are 
often quoted in frequentist statistical interpretations of data. Bayesian statistics turn these 
quantities into probability distributions on the input parameters of the model, requiring 
the introduction of the infamous prior probability distribution. The probability distribu- 
tion of some parameter after confrontation with data is called the posterior probability 
distribution. A global fit of some model to data often consists of finding the variation of 
the figure of merit with the model parameters. The best-fit set of model parameters is 
sometimes quoted, with the amount of parameter space contained within some expected 
amount of statistical variation of data. More complete analyses map out the figure of merit 
on the parameter space, and Bayesian analyses then make probabilistic inferences based 
upon the map. 

7.1 Algorithms for multi-dimensional fits 

Even if one restricts the MSSM to some lower number-of-parameters form such as mSUGRA, 
the parameter space is still of considerable dimensionality: 4 for a given sign of p (rug, Mi / 2 , Aq 
and tan/3). Also, if one wants to perform global fits of the model to data, one should in- 
clude variations of the relevant Standard Model input parameters, mt is proportional to 
the largest parameter in the model, the top Yukawa coupling, and for high tan j3 the bot- 
tom Yukawa coupling, proportional to the bottom quark mass, can change the predicted 
values of observables. Variations of a s within its empirical uncertainties can also have a 
large effect on squark and gluino masses through the RGE evolution, since it is the largest 
gauge coupling. If one is including precision electroweak observables in the fit, including 
uncertainties on the fine structure constant a becomes essential. Thus, in mSUGRA one 
has an eight-dimensional relevant parameter space. Scans in such a space are impractical, 
since the required number of points is exponential in the number of parameters. If one 
required a resolution of 25 points for each parameter, 1.5xlO n points would be required 
in total. To make matters worse, there are often sharp features in the x 2 distribution that 
would render such a low resolution insufficient. Such a large number of points cannot be 
calculated in a reasonable amount of CPU time, even given recent advances in computer 
technology. If one has access to a computer farm, calculating a few million points is feasible 
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within a few days, for example (unless one wants to simulate event generation, which would 
take much longer). There is therefore a need for more sophisticated scanning algorithms 
that can reduce the required number of scanned points for parameter spaces of more than 
three dimensions. 

The software tool MINUIT [§3| is a well-tried function minimiser It calculates deriva- 
tives of the figure of merit with respect to input parameters and performs hill-climbing 
algorithms to try to find the best-fit point. It then determines the error matrix from a 
matrix of second derivatives of x 2 . This error matrix contains information about the la 
standard deviations of the parameters in the Gaussian approximation (where ax 2 is as- 
sumed to be parabolic around the best-fit point) and including correlations. For cases 
where the Gaussian approximation is a bad approximation, another internal MINUIT algo- 
rithm can be used for determining errors including non-linearities, but can be very time 
consuming depending upon the amount of non-linearity. Algorithms that use derivatives 
can be problematic when the surface that they are minimising are rough. In the SUSY 
fitting case, the original SUSY spectrum is obtained by an iterative process up to some 
numerical accuracy, which then feeds into the rest of the figure of merit calculation, pro- 
viding small discontinuities in the surface. A typical numerical fractional accuracy in this 
stage of the calculation might be 1CP 3 . While a fractional accuracy of 10 -5 is feasible, it 
requires much more CPU time per scanned point, and is actually unattainable in certain 
"difficult" regions of parameter space such as the focus-point region. MINUIT also finds 
parameter degeneracies problematic, where the figure of merit does not change much along 
some curve in parameter space. Despite these short-comings, MINUIT has been used to 
perform global fits of mSUGRA to global data successfully fl54| . 

MCMC methods are commonly used in cosmological [5(| and other contexts, and 
recently there has been a realisation that they are very useful to the SUSY high-dimensional 
scanning problem. MCMCs scan more often where the fit is good and the figure of merit 
is high and less often in the tails of distributions. In fact, the density of scanning is 
proportional to the figure of merit. MCMC methods have a high CPU overhead, meaning 
that they are not the most efficient tool for one or two dimensional problems. But the 
required number of points goes roughly linearly with the number of dimensions rather 
than exponentially, and so they are very useful for our higher-dimensional mSUGRA fitting 
problem. In this context, a Markov chain consists of a long list, or "chain" of points and 
their associated likelihoods. Statistical inference can be made by binning these points in 
terms of some quantity of interest. The simplest implementation of MCMC is called the 
Metropolis algorithm [57]. In the Metropolis algorithm, for the first point in the chain, a 
point xo is picked at random in parameter space and its posterior density calculated, pq. A 
potential next point x\ is picked in the vicinity of the previous point, again at random. If 
Pi > POi the new point is accepted. Otherwise, the new point is accepted with probability 
Pi/po- If the new point is not accepted, the previous point is added again on to the chain. 
This algorithm is repeated many times, until it has explored all of the relevant parameter 
space. There are many choices of how to pick a potential next point "in the vicinity" of the 
current one, and some trial and error is usually involved in setting the length scales involved. 
Usually, a Gaussian function is used to randomly choose the distance of the new point away 
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from the current one, but formally, any well behaved function would work in the limit of 
an infinite number of MCMC steps provided it has no true zeroes. For efficient scanning, 
the length scale should be of order the length scale of the likelihood variation. If it is much 
larger, hardly any new points will be accepted and the efficiency will be too low. If it is 
much smaller, many new points will be required to explore all of the good-fit parameter 
space. In order to verify that the algorithm has indeed explored the parameter space 
properly, it is good practice to run several statistically independent chains concurrently. 
One can then compare the results in the different chains statistically to see how similar they 
The Metropolis algorithm does not rely on derivatives and is therefore immune 



arc 



to serious problems caused by roughness from numerical error. It can easily be used to 
interpret data in a Bayesian form or in a frequentist form. For the Bayesian inference, one 
plots the quantities in question (say, squark mass vs gluino mass) in bins. The marginal 
posterior probability distribution in terms of these parameters is then proportional to the 
number of points in the chain that land in each bin. Marginal refers to the fact that all 
other parameters have been integrated over. In order to interpret the chain in a frequentist 
fashion, one plots the profile likelihood: the likelihood of the maximum likelihood point 
that lands in each bin |)9| . Such a procedure, provided a sufficient number of samples to 
get near the maximum for each bin has been obtained, is equivalent to minimising x 2 h" 1 
each bin. Confidence limits can be found in the parameter plane in question by plotting 
iso-A^ 2 contours, where Ax 2 is x 2 assigned to each bin minus the x 2 °f the minimum bin. 
MCMC methods thus provide full maps of the figure of merit across parameter space or 
other scalar quantities that one is interested in. A package SuperBayeS (fortran77, 
f ortran90 and C++) is available for performing global fits to SUSY models using MCMC 



and SOFTSUSY 0, DarkSUSY g§ and FEYNHIGGS @. The MCMC routines were adapted 



from cosmomc ]56j, as well as some of the plotting routines. The program SFITTER |60[ is 
currently being developed which will fit SUSY models to collider data on sparticle masses 
using MCMC methods. 

A problem that is not addressed by either MINUIT or by the Metropolis algorithm 
is that of well-separated x 2 minima. MINUIT only finds a local minima. In principle, 
the Metropolis algorithm may find all local x 2 minima in the limit of infinite number of 
samples. In practice however, if the local minima are small and require small length scales 
for suggesting proposed points, and the distance in parameter space between them is large, 
the chance to "hop" from one local minimum to the other may be tiny and require an 
unfeasibly large number of samples. A "tweak" to the Metropolis algorithm exists which 
can solve this problem and is called bank sampling f6l| . In bank sampling, one performs 
a two-step process. In the first step, many different Metropolis chains are started and run 
for a small number of steps, but numerous enough to find points somewhere near local 
likelihood maxima. These points then form the "bank" or "cache" of points used in a new 
modified Metropolis algorithm. On each MCMC step, there is a small probability that 
the chain will propose a point in the vicinity of one of the bank points. If the new point 
is added successfully to the chain, the chain "teleports" to the other local maxima. In 
this way, the relevant local maxima all appear in the fit results, correctly normalised with 
respect to each other. 
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If only the global likelihood maximum, or equivalently, the global x 2 minimum, is 
desired, a different modified Metropolis algorithm called simulated annealing can be used. 
Simulated annealed can be used to find a point near a global x 2 minimum when several 
local ones exist. In simulated annealing, it is imagined that the x 2 surface is some potential 
energy surface upon which a particle moves. A finite temperature is set, which increases 
the length scale of the proposal step (or, in the analogy, the average distance the particle 
moves) . The temperature is very large at the start of the algorithm and gradually decreases 
to one thereafter. The chance of acceptance of a worse-fit point is also fixed to be higher 
with increased temperature T, being set to e~^ x2 / T , where Ax 2 is the x 2 difference between 
the current point and the proposed worst-fit one. In the early stages, the algorithm is more 
likely to traverse bad-fit regions and not be trapped in local minima. The computer code 
FITTINO [62] can fit a 24-parameter simplified weak-scale MSSM to assumed cross-section 
and mass from SUSY signal collider data. Tree-level values of observables and subsets of 
SUSY parameters are used to obtain start values for the x 2 ~fit- Simulated annealing is then 
performed in order to find a better approximation to the global x 2 minimum. Using these 
parameter values, MINUIT is performed in order to minimise x 2 more precisely. In order to 
investigate the uncertainties in the fit, a series of fits for many imagined experimental data 
are performed in FITTINO, with data smeared around their nominal values, and the global 
X 2 minimum is found in each case. 

In frequentist statistics, hypothesis testing often reduces to finding the minimum x 2 
of different models. However, in Bayesian statistics, one wishes to calculate the evidence 
ratio: the ratio of volumes under the posterior probability surfaces, a quantity that can 
be very computationally intensive to calculate. Bank sampling provides a method for the 
rough computation of the evidence ratio, by having bank points within each of the separate 
models. After the MCMC has run, the ratio of points in each model is an estimate of the 
evidence ratio. Such an estimate may not be very accurate, particularly where the evidence 
ratio is much larger or smaller than one. In such cases, one can artificially multiply one of 
the model's likelihoods by a factor which will bring the resulting evidence ratio closer to 
one. The normalisation can be un-done, with the result that the ratio can be computed 
with smaller statistical uncertainty from the likelihood re-scaling. The disadvantage of 
bank sampling for Bayesian evidence evaluation is that only ratios of the evidence can be 
determined, not the evidence value on its own. 

An algorithm which solves this problem as well as the well-separated likelihood maxima 
problem in a completely different way is the 'MultiNest' technique MultiNest models 
the multi-dimensional likelihood surface with a series of (possibly overlapping) ellipsoids. 
Clustering algorithms are contained within the larger algorithm. They determine when an 
ellipsoid is to be broken up into two different ellipsoids because the initial one does not 
model the underlying distribution well enough. Many live points are chosen, sampled from 
the prior probability distribution. Current live points are described in terms of ellipsoids, 
determined by the covariance matrix of the live points, enlarged by about 20% to take non- 
linearities into account. The live point with the smallest likelihood is replaced with one with 
a higher likelihood re-sampled from the ellipsoids. Thus, the live points gradually home 
in on the likelihood maxima as the algorithm proceeds and the evidence can be calculated 
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from the list of live points and their evidences, as can posterior probability inferences. The 
evidence of a single model can be accurately calculated in this approach, in contrast to the 
case of bank sampling. Thus, as one builds up a list of different models that one is testing 
against some set of data, there is no need to run many different comparisons between the 
different pairs of models: a single computation for each model suffices. The nested approach 
does need to be able to sample efficiently from the prior probability distribution and so will 
not work efficiently in cases where there is no analytic form for the prior. For extremely 
high dimensional cases (say, 10 and above), a MCMC-hybrid nested sampling approach 
may be more efficient than the ellipsoidal approach ||J for multi-modal distributions. 

7.2 Example Global Fits to mSUGRA 

We now display some example results using the various techniques introduced in section |7.1| . 
We pick examples of global fits to mSUGRA in the literature as an example. Typically, the 
data that authors have chosen to fit to include the relic density of dark matter to the relic 
density of neutralinos, M\y, sin 2 0^ , BR(b — > 57), BR(B S — > — m t , m^m^), 

a s (Mz), a(Mz) and direct exclusion limits from colliders. Fig. |]a shows the posterior 
probability distribution in terms of the mo — Mi/ 2 plane for both signs of /i with flat priors 
in mo, M1/2, tan/3 and Aq after an MCMC fit using bank sampling [59|. The probability 



relative to the one in the maximum-posterior bin is shown by the colour, and measured by 
the bar on the right. The most probable region at low values of mo and M1/2 corresponds 
to the stau co-annihilation region, where the lightest stau and lightest neutralino are quasi- 
mass degenerate. The extended probability mass at mo,-/Vf 1 / 2 ~ 0.5 TeV corresponds to 
the ^4° boson resonance at high values of tan j3, where dark matter annihilation proceeds 
efficiently through x?X? ~~ * A ~^ bb. At larger values of mo, we have the focus point 
region, where efficient annihilation into weak gauge boson pairs is possible. In Fig. |]b, the 
same data is interpreted in a frequentist fashion using the profile likelihood technique. This 
technique picks out the best-fit points, rather than averaging over all points in the unseen 
dimensions. Figs. [|a, b differ where there are significant volume effects, that is where the 
volume of points in the unseen dimensions enhances or diminishes the Bayesian fit. The 
fact that the frequentist interpretation differs from the Bayesian one can be seen as a signal 
that more data is required for the fit; indeed we should not be surprised since a complex 
model with eight free parameters has been fit with some fairly indirect data. Similar fits to 
the Bayesian ones above were performed using the Metropolis MCMC algorithm, resulting 
in quite similar posterior probability densities for the particle physics properties, despite 



some differences in the indirect constraints used [64]. In addition, ref. [64] constrains dark 
matter detection cross-sections. We show the posterior probability distribution function of 
the spin-independent direct dark matter detection cross-section in Fig. |]c for flat priors in 
mo, M1/2, tan/3 and Aq, for fi > 0. The most constraining direct detection experiment, 
XENON-10, can be seen to cover some of the favoured region already, assuming that the 
flux of dark matter passing through XENON-10 is the same as the galactic average. In 
Fig. 0d, we see the results of a more traditional frequentist x 2 mSUGRA fit using MINUIT 



in terms of the lightest CP-even Higgs mass of the MSSM [54]. For each value of the 
lightest Higgs mass, a x 2 minimisation was performed against all of the other mSUGRA 



- 19 - 




Figure 4: Global mSUGRA fits in the mo — M1/2 plane: (a) shows the Bayesian posterior proba- 
bility distribution j5J|, (b) shows the frequentist interpretation in the same plane j5jj], (c) displays 
the direct spin independent detection cross section posterior probability distribution function ver- 
sus mass of the lightest neutralino for \i > along with some 95% C.L. exclusion contours from 
direct detection experiments |64| . Inner and outer contours show the 68% and 95% confidence level 
regions respectively, (d) shows the A\ 2 of the lightest CP-even Higgs mass from all constraints 
except for the direct LEP2 Higgs mass constraint |M. 



parameters. Many additional electroweak and b observables were included in the x 2 °f this 
fit, although for comparative purposes the LEP2 direct bound on the Higgs mass was left 
out. This bound is plotted as the yellow excluded region in Fig. ||d. It can be seen that 
the global x 2 minimum occurs at ~110 GeV, below the direct 95% C.L. lower bound of 



114.4 GeV. The authors of Ref. |54j use the x 2 curve to infer that m h = 110jlf ± 3 GeV, 



where the first uncertainty is statistical and the second uncertainty is theoretical. 
7.3 Data archival 

The samples from MCMC fits take some effort and CPU time to obtain. In principle, the 
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mSUGRA fits could be useful to other physicists, who wish to make their own inferences 
about observables. If a new calculation of SUSY contributions to some observable were 
to come on-line, statistical inference could be made by simply obtaining some indepen- 
dent samples from the MCMC chains and calculating the new observable for each. A few 
thousand points might suffice in terms of statistics. Then, the totality of current empirical 
knowledge about SUSY corrections to the observable is obtained without a need for com- 
plicated multi-dimensional fitting procedures. Aside from that, other physicists might be 
interested in using the chains for their own scans over the points in parameter space that 
are compatible with current data. For this reason, the authors of Ref. have formed the 
KISMET web-site, which contains links to text files of the chain data. The weight of each 
point (the number of times it was visited in the MCMC procedure), along with the values 
of input parameters, resulting indirect observables, sparticle masses and the likelihood, are 
listed in the files. Also, 10 000 independent samplings from the chains in SLHA format are 
available from the web-site. 



The SuperBayeS [52] data are now available to some extent on-line: one can fill out 
a web-form in order to automatically receive plots of the posterior probability density in 
dimensions specified by the user (65]]. These dimensions can be specified to be observables, 
input parameters or sparticle masses. 



8. Summary 

Currently, little is empirically known about supersymmetry except for a few indirect data. 
This tends to lead to under-constrained SUSY models and consequent degeneracies in 
global fits. However, if the LHC provides some signals that are compatible with SUSY, 
aside from being an extremely exciting discovery beyond the Standard Model, hypothesis 
tests against alternative models and even between different classes of SUSY models will be 
desirable. Ideally, constraints upon the SUSY Lagrangian would be derived with the help 
of high energy e + e~ linear collider experiments. A well defined theoretical framework is 
needed when higher order corrections are included in trying to reconstruct a fundamental 
SUSY theory and its breaking mechanism. For this purpose, the Supersymmetry Parameter 
Analysis (SPA) |66|] scheme provides a consistent set of conventions and input parameters, 
as well as a repository for programs which connect parameters in different schemes and 
relate the Lagrangian parameters to quantities that may be more directly extracted from 
physical observables such as masses, mixings, decay widths and production cross sections 
for supersymmetric particles. 

There is a somewhat bewildering proliferation of computer program tools for SUSY 
calculations and phenomenology in the literature. This proliferation is a useful develop- 
ment, and reflects the interest of the high energy physics community in supersymmetry. 
Even more useful is the collusion, collaboration and organisation between the different 
programs, to allow results from one to be fed into another program and interpreted auto- 
matically. The SUSY Les Houches Accord is a good example of such practice, and it has 
now become essential for any relevant computer tool to use it so that it can communicate 
with the other tools. The most popular computer languages that the tools are written in 
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are still various versions of f ortran and C++. Following the move of most new high energy 
physics experiments to C++, there is a tendency for new event generators to be written in 
C++ rather than f ortran. In fact the precise language a SUSY tool is written in is be- 
coming less important with the advent of the communication accords, which are in ASCII 
format. Many of the more sophisticated matrix element or event generators use their own 
encoding of a Lagrangian to enable the user to define new models. This approach is of 
obvious use and generality, but packages have many different definitions of the Lagrangian. 
Perhaps there is a need for yet another accord, so that the same model can easily be fed 
in to different tools without the need for the user to translate the Lagrangian between the 
various different conventions. 

Finally, we end with a brief summary of the phenomenological SUSY tools that are 
covered in this review. Since some programs have several different functions, we summarise 
them all together in Table |2[ although the programs are loosely grouped according to their 
functionality. C indicates that the tool includes or uses a method of encoding a Lagrangian 
in order to define extended or new models. In the table, 'Spectrum' indicates that the 
tool includes an SUSY spectrum calculator, vr indicates that RGEs include an option 
for including right-handed neutrinos (and therefore neutrino mass models), RPV indicates 
that the tool can handle R— parity violation, NMSSM that it can calculate in the Next- 
to-Minimal Super symmetric Standard Model, CPV that the tool can take into account 
complex phases in the SUSY sector and FV that the tool includes some non-minimal flavour 
violating effects. 'Decays' indicates that the tool automatically calculates the branching 
ratios of SUSY or SUSY Higgs decays in the MSSM or extensions. Tools which have a 
positive entry under 'Decay spin' include angular correlation effects from sparticle spins 
when simulating decays down cascade decay chains. 'ME' indicates a matrix element 
generator: the code can simulate scattering for 2^ N hard particles, where N > 2. 'Initial 
pol' shows that polarisations of the colliding particles can be taken into account: usually in 
e + e~ collisions, but sometimes also in 77 or ej collisions. A tick under the osusy heading 
means that the code has an easy user interface for calculating total or differential cross- 
sections for the production of (sometimes specified) sparticles and/or SUSY Higgs. e + e~ 
and pp indicates that the initial colliding particles can be leptonic or hadronic, respectively. 
'Events' mean that individual events are simulated, PS/Had that the program can perform 
parton showering and/or hadronisation of partons. A tick under the QjjMh 2 header means 
that the relic density of dark matter can be calculated, cr p SI SD that an estimate of dark 
matter direct detection is included and 'Ind. DM' that some indirect dark matter detection 
fluxes are provided. For 6-observables, b — * sj, B —* tv, tt and B s — > p + p~ indicates that 
there is a calculation of the relevant branching ratio including some SUSY effects. A 
positive entry for Ao- means that the program calculates the isospin asymmetry in B 
decays, whereas an entry under AMb s that the SUSY contributions to B® d — B® d mixing 
are calculated. A {g — 2)^ entry indicates that a SUSY contribution to the anomalous 
magnetic moment of the muon can be easily extracted from the tool, whereas EW means 
that some electroweak observables are provided: usually Ap and M\y. A tick under edm 
means that electric dipole moments can be calculated, whereas 'Fits' indicates a fitting tool 
that can fit either collider observables and/or indirect constraints such as EW observables 
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and dark matter relic densities. An entry under 'Web form' gives the reference including a 
link to a web-form where results from the program can be automatically obtained by filling 
in a form on the world-wide web. Finally the 'code' column indicates that the package 
can output computer code, which can then be compiled into a numerical program in order 
to evaluate observables. Prospective users are warned that multi-functionality does not 
necessarily mean a more accurate calculation and indeed in some cases, the converse will 
apply. It is hoped that Table || will help point prospective new users towards the SUSY 
tool(s) that they require. 
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[17] 


V 


SUSYGEN [37] 






V 




V 








V 


V 




V 




V 
































Wliizard [35] 


V 


















V 




V 


V 


V 
































Herwig++ [16] 


V 




V 








V 


V 






V 


V 


s/ 


V 


V 






























ISATDDLS [9] 




V 










V 






V 




V 


V 




V 


V 


V 




V 


V 






V 


V 








[14] 




PYTHIA [18] 








V 






V 










V 


V 




V 






























SHERPA [33] 


V 








V 


V 




V 


V 


V 


V 


V 


</ 




V 






























_|ter«USY [48JI 

micrOMEGAs [49] 
PRDSPIND [40] 


D 


D 


J 


D 

V 


D 

V 


V 






D 


V 


D 

V 


D 

V 






V 


V 

•J 


V 
V 


V 
V 


V 




D 


J 


V 
V 


V 


V 




[48] 
[14] 




Superlso [51] 
SusyBSG [50] 




V 










V 




V 






V 












FITTIND [62] 
KISMET [59] 
SuperBayeS [52] 




D 
D 
D 


D 
D 








D 


D 


D 






D 
D 


D 


D 


J 
V 


J 








J 

D 


J 

D 




V 
V 


[59] 
[65] 





Table 2: Summary of functionality of current, publicly available, supported SUSY tools. A y/ indicates that there is some support for the feature 
in question, but makes no claims about the accuracy of the calculation. D indicates that the one of the other packages in the table is included in 
the distribution in order to provide the relevant functionality. See section 8 for a description of the various features. 



